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Abstract 

This paper describes a method to filter oscillatory transients from mea- 
surements of a time series which were at least an order of magnitude larger 
than the signal to be measured. Based on a Kalman filter, it has an opti- 
mality property and a natural scaling parameter that allows to tune it to high 
resolution or low noise. 



1 How transients can spoil measurements 

In her Matura thesis [[5]|, Carmen Miiller tried to measure the opening shock of 
small parachutes in the Ackeret wind tunnel of the Federal Institute of Technology 
in Ziirich (ETHZ). A specially constructed mounting tube ejected the parachute 
using a piston powered by a small black powder charge into the air stream. While 
this accurately reflected the way parachutes are usually ejected from small rock- 
ets, it also induced heavy oscillations in the measurement equipment and probably 
also excited waves within the wind tunnel. As a result, the data was almost unus- 
able. The large amplitude transients were at least one order of magnitude larger 
than the effects to be measured (figure [TJ. 

Upon inspection of the data, two oscillatory components with frequencies of 
about 40Hz and 12.2Hz were found. Assuming that the measured signal is a 
superposition of the actual force to be measured f(t) and a transient of period 
T, convolution with a function of constant value T^^ on the interval [— T, 0], or 
equivalently a moving average over the samples that fit in an interval of length 
T, eliminated the transient. In this way, the two dominant transients could be 
removed, and an approximation to the opening shock force was obtained. The 
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Figure 1 : Wind tunnel measurements of opening shock dominated by transients 
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Figure 2: Opening shock measurements filtered with moving averages of different 
lenghts. The red curve obviously does not accurately represent the forces during 
parachute opening. 
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price to pay was of course that the signal / was also modified by the convolution. 
The modifications were large enough that it became unreasonable to draw any 
quantitative conclusions from the measurements (figure [2]). 

Using an exponential weighting function, it was possible to slightly improve 
the filtering, taking into account the damping of the transients. Since the carrier 
of the convolution of the intervals of length 24ms and 83ms corresponding to the 
frequencies of the transients is as an interval of at least 107ms, the convolution 
reduced the effective resolution of the measurements to less than 10 samples per 
second, which is of the same order of magnitude as the force changes during the 
opening shock. 

While this approach was successfuU und fully adequate for a Matura thesis, the 
question was still open whether a better filtering method could reveal more details 
about the opening shock. Ideally, such a filter would have a scaling parameter that 
allows to tune it to different resolutions, where of course fine grained resolution 
comes at the price of additional noise and incomplete filtering of the observed 
transients. Below we show that a Kalman filter is ideally suited for this situation, 
and that the system errors represent a set of parameters that allow this kind of 
tuning of the filter. 

2 How to get rid of transients 

In the construction below, we will not develop the Kalman filter equations, but 
only indicate the system equations from which the filter can be derived. Our 
notation closely follows standard texts on the discrete Kalman filter like [|3l or Pll. 

2.1 The harmonic oscillator 

The measurement system in the wind tunnel works as follows. The forces on the 
parachute and the ejection tube deflect the mounting hardware by a tiny fraction 
measured by piezoelectric sensors. So the mass m of the mounting hardware is 
subject to two forces, the forces we want to measure and the elastic forces, and 
we measure its position x as a function of time. If D is Hook's constant for the 
elastic forces, x obeys the differential equation 

mx + ax + Dx = f, (1) 

where / is the force we want to measure. In a pseudostatic situation, we can 
assume that x = x = 0, and so the force is proportional to the deflection value we 
measure. In the dynamic situation, we have to infer f{t) from x{t) as a function 
of time. 
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A naive approach would use the measurements of x and compute approxima- 
tions of X and x as finite difference quotients. Then an approximation for / can be 
computed from ([!]). However, numerical differentiation tends to be noisy, and the 
second derivative will probably be almost useless. The problem gets even more 
serious in light of the fact that in real applications, we need to filter several times, 
so we need a method to estimate / that does so with a controlled level of noise. 

This is the classical setting for a discrete Kalman filter. At any point in time we 
would like the best estimate of the current position and speed of the mass m and 
the force / based on measurements of x at discrete points nAt in time. We use 
Xn = {x{nAt), x{nAt), f(nAt)) as the state vector of a linear system. The time 
evolution is given by the differential equation ([T]), which can be written, using the 
familiar substitution v{t) = x{t) in vector form as 
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Since we don't know anything about /, we assume for the purposes of the Kalman 
filter, that / is constant, i. e. / = 0. Using the state vector s = the 
differential equation becomes 

—s = As 
dt 

where 

A = 

Time evolution for a short time interval At is given by the matrix exponential 

if = e^^*. 

The measurement equations are also very simple. Since we only measure the 
position, we can use the measurement matrix 

H = {1 0) . 

To completely define the Kalman filter, we have to specify the system and 
measurement errors. The measurement error is a well known characteristic of the 
measurement equipment. The system errors for x and v have to be kept very small, 
as we believe the harmonic oscillator equation is followed quite accurately. In 
contrast, the evolution equation for the force is only a very rough approximation. 
We attribute almost all deviations of the real system from the model to external 
forces. Thus the flucutations due to measurement errors are the only sources of 
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Frequency 



Figure 3: Frequency response of Kalman filter based transient filter for various 
values of aj 

errors for x and v, while for / we have to accept errors as large as the slope of / 
times the time step. If cr^. is the position measurement variance and cr| the vari- 
ance of the force, then a reasonable approximation of the system error covariance 
matrix Q and the measurement error covariance R is 

/al 0\ 
Q = crl , R= {al) . 
\0 ajj 

In the example of the wind tunnel measurements of the opening shock, the 
measurement error was about O.OOIN^, while the force error was on the order of 

Using a large value of a'j causes the Kalman filter to "trust" the measurements 
and agressively attribute deviations from the evolution given by the differential 
equation to changing external forces. Of course, / will be much more noisy in this 
case, since measurement noise is converted to a large degree to noise in /. If aj 
is small, then the filter tries not to modify / too quickly, and attributes deviations 
from the differential equation to a larger extent to the measurement errors. The 
resulting / will be much smoother, but steep changes in / will be filtered away. 
Thus cr? is a natural scaling parameter we asked for in the introduction. 



2.2 Frequency response 



Figure [3] shows the numerically computed frequency response for various values 
of CTj. As expected, the notch at 50Hz is very pronounced. For large enough 
values of cxj, there is hardly any loss for high frequencies. For small values of a'j, 
however, the loss is much larger, rather similar to the sinx/x characteristic of the 
moving average filter. 



2.3 Implementation considerations 
2.3.1 Building the filter from on 

None of the coefficients m, a and D of the differential equation are known. Only 
the frequency and the damping of the transients can be measured. The character- 
istic equation of the harmonic oscillator equation 

uj^ + —uj + — = (3) 
m m 

has two complex conjugate roots u and u. To simplify the notation, we write 
a = ^ and b = —. We conclude that 



a . , La? 



Re u; = - and Im = \ b — 

2 V 4 

Given uj, the differential equation can be obtained from 

a = 2Reu; and 6 = (Imcj)^ + 



2 «^ 



This allows to easily construct the Kalman filter for the oscillator from frequencies 
obtained e. b. by spectral analysis. 



2.3.2 Precomputing the Kalman gain matrix 

It is well known that when all the matrices defining the Kalman filter are constant, 
then the Kalman gain matrix converges. For given values of the system and mea- 
surement errors, the Kalman gain matrix can thus be computed beforehand and 
"hard wired" into the filter. 



2.4 Coupled oscillators 

If a system exhibits several oscillations, we can still assume that the system is 
linear, and that each oscillation corresponds to a harmonic oscillator with that 
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Figure 4: Harmonic oscillator excited by a pulse of 0.1s duration at time t = 1. 

frequency. Such a system can be described by a higher dimensional state vector, 
and it is possible to construct a Kalman filter for it. However, it would become 
rather large, and the following simpler approach is preferable. 

Each one of the coupled oscillators is accurately modelled by a harmonic os- 
cillator on which two external forces act, the force / and the combined forces 
excerted by the other oscillators. The latter have a well known frequency, so fil- 
tering the signal in turn for all the frequencies of transients identified using the 
method described above removes all the transients from the signal /. 

3 Success stories 

The method described in the previous section was first applied to synthetic data to 
verify that the filter is indeed capable of retrieving the synthetic signal. Then the 
data was applied to the opening shock data that originally motivated this research, 
and to a related problem involving thrust curves of rocket motors. 

3.1 Test data 

To verify the properties of the Kalman filter based deconvolution of transients, a 
solution for the harmonic oscillator equation 

X + 0.1- X + 1000- X = f{t) 
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Figure 5: Signal from figure |4] filtered by a Kalman filter for increasing aj, com- 
pared to the original excitation (blue curve). 

was computed numerically for the force function 

Jo t < 1 Vt > 1.1 

1 10000 tG [1,1.1] 

and some gaussian noise with a = 0.1 added to simulate the measurement errors. 
The resulting measurement data is shown in figure|4j Figure [5] shows the output of 
a Kalman filter based deconvolution with a'j = 10'', k = 3,4:, 5, 6. As expected, 
noise increases with increasing a'j. However, a simple moving average over just a 
few samples can remove large parts of the noise without distorting the signal very 
much. 



3.2 Opening shock of parachutes 

As a practical test case, the filtering method was applied to the opening shock data. 
A chain of three kalman filters was constructed using the following algorithm. 

1 . Take the raw as the current signal to filter. 

2. Find the dominant frequency to in the current signal. 

3. Construct a Kalman filter for frequency to. 
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Figure 6: Opening shock measurements from figure [T] filtered with three Kalman 
filters for different values of the force error variance a'j. 

4. Filter the current signal, increase a'j until a new dominant oscillation ap- 
pears, in which case you continue at step 2, or the signal becomes implau- 
sible. 

5. Reduce until you again get a physically plausible signal. 

During this process, in addition to the frequencies already found in f5\, an addi- 
tional oscillation at 85.7Hz was found. The result of the filter for different values 
of aj is shown in figure [6| The curves following most closely the physical force 
signal were obtained with aj = lOOOOON^ and aj = 200000N^ the second again 
shows some probably unphysical transients. 

The reconstruction of the signal using the Kalman filter at least allows now to 
draw some quantitative conclusions from the opening shock data. E. g. it allowed 
to confirm the empirical rule used in some circles that the opening shock is about 
twice as strong as the drag of the fully inflated parachute. 

3.3 Thrust curve of high thrust rocket motors 

Aerotech Consumer Aerospace produces rocket motors for hobby rockets and 
small research rockets. For each motor, thrust curves are available from various 
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Figure 7: Thrust curves of 38mm Warp-9 propellant rocket motors as they appear 
on the instruction leaflet accompanying the reload kits. H999 top curve in top 
diagram, 11299 in bottom diagram. 



Flight 2, 29. 1. 2008, H669P 
boost phase 



Figure 8: Acceleration profile of a rocket flight on a H669 motor, no oscillations 



10 




O.OOS 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 O.OS 0.0S5 0.07 0.075 0.03 0.085 0.09 



Figure 9: Original thrust curve (blue) and version filtered for transients using a 
Kalman filter (green) 

sources on the Internet, and also from Aerotech. The more recent motors using 
on the Warp-9 propellant have extremely high thrust for a very short time interval. 
The motor H999 has average thrust of 999N during only 0.3s. The thrust curve 
shown in figure |7] bottom (from [1]) shows an initial oscillation which cannot 
be explained by any known mechanism in solid propellant motors. As a matter 
of fact, measuring the acceleration during an actual rocket flight with this mo- 
tor, as was done with the high precision altimeter developed during the rocket- 
ware project (http://rocketware.othello.ch), shows no such oscillation (see figure 
[8]). Also the only slightly larger Motor 11299 (figure |7] bottom, from [2J) has a 
completely flat thrust curve. So most probably the oscillation is an artefact of 
the measuring equipment used for these curves. Of course, measuring thrust of a 
rocket motor is similar to measuring aerodynamic forces in a wind tunnel. Thrust 
acts on the motor, which is displaced slightly, the displacement is proportional to 
the force. But without sufficient damping, the test stand holding the motor may 
oscillate adding an oscillatory transient. 

To corroborate this hypothesis, the Kalman filter approach was applied to a 
manually digitized version of the thrust curve. Figure |9] shows the filtered version 
which is obviously much closer to the thrust curve measured by the rocketware ac- 
celerometer. The remaining noise was probably introduced when manually read- 
ing thrust values from the printed thrust curve. 

It turns out that the thrust curves were in fact measured on different test stands. 
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Thrust curves for 11299 were obtained on Aerotechs test stand, the test stand of the 
Tripoli Rocketry Association was used for motors that show strong oscillations. 

4 Conclusion 

The results of actual measurement data filtered using the Kalman filter approach 
show that oscillatory transients of at last the magnitude of the signal can be ef- 
fectively filtered to reveal the signal. The measurements must be precise enough 
to render the full dynamic range, otherwise the noise introduced by the filter will 
quickly become larger than the signal. If this condition is satisfied, however, the 
method can be expected to be close to optimal, as the Kalman filter has such a 
property. If better models for / are available, it is possible that the method could 
be refined to yield even better results. 
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